function obj_eval = NLI_obj_funct(param,fp,data,col_long,NLI_data_mom,estimate_NLI_process)

% Baseline model for non-labor income process: two-step logit

% 1) logit parameters
mu1 = param(1); %intercept
mu2 = exp(param(2)); % mom_age: has to be positive
mu3 = param(3)-1; % mom_age squared
mu4 = exp(param(4)); % mom_educ: has to be positive
mu5 = exp(param(5)); % dad_age: has to be positive
mu6 = param(6)-1; % dad_age squared
mu7 = exp(param(7)); % dad_educ: has to be positive

% 2) step two parameters
mu8 = param(8); %intercept
mu9 = exp(param(9)); % mom_age: has to be positive
mu10 = param(10)-1; % mom_age squared
mu11 = param(11)-1; % mom_educ
mu12 = exp(param(12)); % dad_age: has to be positive
mu13 = param(13)-1; % dad_age squared
mu14 = param(14)-1; % dad_educ
sigma_NLI = exp(param(15)); % shock std: has to be positive

vp = [mu1;mu2;mu3;mu4;mu5;mu6;mu7;mu8;mu9;mu10;mu11;mu12;mu13;mu14;sigma_NLI]; % save vector of parameters

% simulate a matrix of non-labor incomes at the optimal parameters
[NLI_sim_data, NLI_sim_data_3d] = NLI_sim_funct(vp,fp,col_long,data);

% Note: censoring and truncating has already occurred in NLI_sim_funct
NLI_sim_mom = NLI_moment_funct(NLI_sim_data,col_long); % calculate simulated moments

diff = NLI_sim_mom - NLI_data_mom;
diff = diff(fp.NLI_moments_with_variation);
diff(isnan(diff)) = 0;
goal_fct = diff'*fp.NLI_wt_matrix*diff;

% if fp.estimate_NLI_process == 1
if estimate_NLI_process == 1
    obj_eval = goal_fct;
else
    obj_eval.NLI_param = vp; % transfer optimal "transformed" parameters
    obj_eval.NLI_param_untrans = param; % transfer optimal "UNTRANSFORMED" parameters (will be needed in bootstrap_SE.m)
    obj_eval.NLI_sim_data = NLI_sim_data; % large MATRIX
    obj_eval.NLI_sim_data_3d = NLI_sim_data_3d; % large MATRIX
    obj_eval.NLI_sim_mom = NLI_sim_mom; % all simulated moments
    obj_eval.NLI_data_mom = NLI_data_mom; % all data moments
    obj_eval.obj_value = goal_fct;
    temp = zeros(length(NLI_data_mom),1);
    for j = 1: length(NLI_data_mom)
        temp(j) = (NLI_sim_mom(j) - NLI_data_mom(j))^2 * fp.NLI_wt_matrix(j,j);
    end
    obj_eval.NLI_weight_perc = 100 * temp/sum(temp); % gives the percentage of obj_eval caused by moment j
    obj_eval.NLI_wt_matrix = fp.NLI_wt_matrix;
    obj_eval.NLI_moments_with_variation = fp.NLI_moments_with_variation;
end